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Abstract 



Quantum-Mechanical methods that are both computationally fast 
and accurate are not yet available for electronic excitations having 
charge transfer character. In this work, we present a significant step 
forward towards this goal for those charge transfer excitations that 
take place between non-covalently bound molecules. In particular, 
we present a method that scales linearly with the number of non- 
covalently bound molecules in the system and is based on a two- 
pronged approach: The molecular electronic structure of broken-symmetry 
charge-localized states is obtained with the Frozen Density Embedding 
formulation of subsystem Density-Functional Theory; subsequently, in 
a post-SCF calculation, the full-electron Hamiltonian and overlap ma- 
trix elements among the charge-localized states are evaluated with an 
algorithm which takes full advantage of the subsystem DFT density 
partitioning technique. The method is benchmarked against Coupled- 
Cluster calculations and achieves chemical accuracy for the systems 
considered for intermolecular separations ranging from hydrogen-bond 
distances to tens of Angstroms. Numerical examples are provided 
for molecular clusters comprised of up to 56 non-covalently bound 
molecules. 
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1 Introduction 



Charge transfer (CT) reactions are ubiquitous in all branches of chemistry, 
and huge efforts have been spent in the past to analyze these processes in 
terms of theoretical models [IH5]. They can be categorized into process in 
which a charge separation in a donor-acceptor system takes place, 

D + A -> D + + A~, (1) 

and processes where the charge has been externally created (either by inject- 
ing or removing an electron), 

D- + A -»■ D + A~, (2) 
D+ + A ->■ D + A + } (3) 

where we define D as the donor and A as the acceptor. For the reactions in 
dlj) and donor and acceptor are defined in terms of donating and receiving 
electrons, while in (j3J) the definitions are in terms of the location of the excess 
positive charge. The three kinds of CT reactions in (JTH3I) are involved in a 
plethora of important processes. For example (p^) may resemble a charge 
splitting event either at a semiconductor interface, between an organic dye 
and a semiconductor, or the charge separation event in the reaction centers 
of photosynthetic systems. Generally, all the charge separation events can be 
thought of in terms of a two-state model comprised of a neutral state and a 
charge separated state. Throughout this work, we call the reactions involving 
only the motion of an excess of charge, as in ([2]) and (j3J), migration CT 
(MCT, hereafter) reactions. MCT reactions are the simplest type of reactions 
involving CT states, and take place in many fundamental processes related 
to charge mobility, transport, and enzyme functionality. In these processes, 
the charge separation event has happened in the past, and the problem shifts 
to the prediction of the kinetics of the excess charge, may that be across a 
polymer [1H6] (organic electronics) , a protein (enzyme functionality [7J[8] and 
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photosynthesis [SIHD]), or DNA [TTHT3] (oxidative damage). 

Very often, the ground and MCT states are close in energy and one easily 
finds the first MCT state as being the first excited state of the system. This 
is an important advantage in theoretical studies of this type of states as 
opposed to other types of excited states for which the excited state search 
can be a daunting task. Another important quality of MCT states is the 
simple physical depiction of ground and excited states: the former has an 
excess charge on the, say, "left" , and the latter features an excess charge on 
the "right". Obviously, this simple depiction becomes more complicated as 
soon as the system under study contains more than two chemical moieties 
capable of accepting/donating an excess charge. These unique qualities make 
MCT states a perfect testbed for new computational methods aiming at the 
accurate prediction of generally all kinds of CT excitations. 

There are two fundamentally different schemes for calculating CT exci- 
tations. The first one involves the construction of Hamiltonian and overlap 
matrices in a basis of charge-localized broken-symmetry states. This basis 
set is often called a diabatic basis set, and the basis functions termed dia- 
batic states. The chosen broken-symmetry states in the set should resemble 
the true ground and the CT excited states of the system. Depending on the 
size of the basis set (or equivalently the number of broken-symmetry states 
included in the calculation) the simultaneous diagonalization of overlap and 
Hamiltonian matrices (i.e. solving the generalized eigenvalue problem) yields 
the sought ground and CT excited states. The latter diagonalization step 
is similar to the diagonalization step in the Configuration Interaction (CI) 
method carried out with valence bond configuration state functions. Even 
though this is a general method, it is particularly suited for the calculation 
of CT excited states. This is because contrary to valence excitations, the 
broken-symmetry basis functions assume a simple charge-localized charac- 
ter and relatively uncomplicated computational criteria can be devised for 
their generation [T4"] - f20] . This is the method that we use in this work and 
is described in detail in the next section for a basis set comprised of two 
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broken-symmetry states. 

Contrary to the first scheme mentioned above, the alternative is not tai- 
lored specifically to CT excitations and it involves obtaining directly the 
adiabatic ground and excited states of the system. As such it is more general 
than the first scheme and is most often adopted in the literature. Due to the 
high computational throughput of modern computers, calculations based on 
the latter scheme are now commonplace. A plethora of quantum chemistry 
methods are devoted to the prediction of excited states' energies, transition 
moments, densities, and wave functions. Despite this, CT excited states have 
always been more challenging than others, because for this kind of excitations 
the electron density changes dramatically compared to the ground state one. 
Density-functional theory (DFT), and specifically its time-dependent exten- 
sion (TD-DFT) through the linear-response formalism for the calculation 
of electronic excitations and transition moments [2TJ has struggled to cor- 
rectly deal with CT excitations when employing approximate density func- 
tional [22]. Pragmatic corrections exist, such as the one originally developed 
by Gritzenko et al. [23j|23], or the one involving the Peach factor [25J. Range 
separated [261428] density functionals have also offered an effective solution 
to the problem by enforcing the correct long-range behavior at the price of 
including yet another parameter to the density functional (the range sepa- 
ration parameter, often denoted by 7). Methodologies to obtain 7 param- 
eters which satisfy certain system-dependent and pro cess- dependent condi- 
tions have been proposed [29|l30j. Novel computational approaches to obtain 
charge-transfer excited states within a DFT formalism are being explored. 
Specifically, a variational formulation of TD-DFT [31] has shown to alleviate 
or to completely cure the CT excitation failures of linear-response TD-DFT 
when the fourth-order relaxed constrained-variational TD-DFT method is 
employed, however, at the expense of higher computational complexity than 
standard linear-response TD-DFT. 

Wave function based methods also have encountered difficulties when ap- 
proaching CT excitations. The low end of wave function methods, configu- 
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ration interaction with singles (CIS) and time-dependent Hartree-Fock (TD- 
HF), are known to grossly overestimate CT excitation energies [32H3D]. More 
balanced methods, such as multi-configuration self-consistent field (MCSCF), 
multi-reference CI (MR-CI), or perturbative corrections to CIS, such as 
CIS(D), as well as ADC(2) [37] are able to deliver good CT excitation en- 
ergies [321ES] and are often taken as benchmark [39]. MCSCF, however, is 
known to fail near avoided crossings (a feature always present in systems 
featuring MCT excitations) unless state averaging is employed [ID]. The 
equation-of-motion Coupled Cluster (EOM-CC) theory [4DJ[4l] is also known 
to fail in reproducing CT excitations with a similar accuracy than ionization 
potentials and electron affinities [I2JH3] due to its non size-extensivity stem- 
ming from the non-exponential from of the EOM-CC ansatz for the excited 
states. Linear-response CC, and specifically an approximation to it including 
only single and double excitations known as CC2 has been very successful 
in predicting CT excitations [HHI7]. A particularly powerful implementa- 
tion of CC2 utilizing the resolution of the identity (RI-CC2) ]48j is routinely 
applied to systems with up to 100 atoms. It is now understood that to ob- 
tain an accurate description of CT excitations the employed method must 
describe the dynamic correlation of the CT excited states similarly to the 
one of the ground state [IPJ. In addition, the computational costs associated 
with all the wave function methods mentioned above (with the exception of 
TD-HF and CIS) are prohibitive for most systems in condensed-phase and 
biosystems of interest. 

In the next section, we will describe in detail how ground and CT excited 
state energies and wave functions can be obtained starting from a basis set 
comprised of two broken-symmetry states. In Section [3] we will introduce the 
Frozen Density Embedding fomulation of subsystem DFT (a linear-scaling, 
full-electron electronic structure method) which we use to determine the elec- 
tronic structure of the broken-symmetry, charge-localized states. Then in sec- 
tion H] we will present the theory behind the calculation of the full-electron 
Hamiltonian and overlap matrix elements among the broken-symmetry states 
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(generated with subsystem DFT) which we have implemented in the Ams- 
terdam Density Functional (ADF) computer software [5U]lo"T]. Section 
[5] is devoted to benchmarking of the method for two selected cases of MCT, 
e.g. hole transfer in water and ethylene dimers at several inter-monomer sepa- 
rations. Section [H] features two pilot calculations of the hole transfer in water 
and ethylene clusters containing up to 56 and 20 molecules, respectively. 

2 Excitation energies from a model of two 
broken-symmetry states 

Consider a system comprised of a charge donor and a charge acceptor moiety 
in the absence of low- lying intermediate bridge states (DA system hereafter). 
Such a system can be well characterized by a two-state formalism. That is, 
it is enough to consider either the adiabatic ground and first excited state or 
a set of two broken-symmetry charge-localized states to capture most of the 
underlying physics. There is a large literature [T5lfT6l[52T[5l] supporting the 
idea that the states with localized excess charges, i.e. one where the charge 
is localized on the donor and, conversely, one where the charge is localized 
on the acceptor, are an appropriate basis for modeling the process of charge 
migration between donor and acceptor. 

A state with localized charges generally is not the ground state of the 
system Hamiltonian. Therefore, in the basis of charge-localized states, the 
Hamiltonian matrix is non-diagonal. In addition, as these states are of 
broken-symmetry character, they are generally non-orthogonal to each other. 
Throughout this work, we will refer to the full-electron wave functions of the 
two charge-localized states as $i and $2- For DA systems, the Hamiltonian 
(H) and overlap matrices (S) in this basis are of dimension 2x2, namely 
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Sl2 




(5) 



The definitions in Eqs. fll]) and ([5]) allow us to obtain the CT excitation 
energy as the energy difference of the orthonormal adiabatc states, which 
can be obtained by solving the generalized eigenvalue problem, 



Hu — E H\2 — ES\2 
H 2 i — ES\2 H 2 2 — E 



(6) 



The above equation yields two energy eigenvalues (say E\ and E 2 ) and their 
difference (AE) is the sought CT excitation. In a closed form, we obtain 



^= x r i r:r (7) 




where V 12 = ^ [H 12 - Su^ 

The problem is then shifted to calculating four matrix elements, i.e. the 
two diagonal Hamiltonian elements (Hu and H22, e.g. the total energies of 
the charge- localized states or diabatic energies), the Hi 2 off-diagonal matrix 
element, and the S12 overlap element. The Vu element introduced above is 
commonly referred to as electronic coupling. 

Several methods and algorithms have been proposed to generate ad-hoc 
charge-localized states and to calculate the above matrix elements both for 
purely wave function methods [rp ^ fTS I H^ llSpT - ISS] as well as DFT methods 
pHIT5l[T71[58|[59HH2]. Regular wave function and DFT methods tend to 
scale non-linearly with the system size. Therefore, important systems of 
biological interest, large organics and hybrid organic-inorganic systems (such 
as dye-sensitized solar cells) are largely out of reach of these methods unless 
truncated model systems or substantial approximations are introduced at the 
electronic structure theory level. 

In this work, we build upon ideas presented in a previous paper [62] 
and we construct the broken-symmetry charge-localized states using a linear- 
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scaling technique called Frozen Density Embedding |63J. The construction 
of the charge-localized, diabatic states is a prerequisite in this formalism for 
the calculation of the CT excitation energies. A brief description of this 
technique follows this section. 



3 The Frozen Density Embedding Formula- 
tion of Subsystem DFT 

Subsystem DFT is a successful alternative to regular Kohn-Sham (KS) DFT 
methods due to its ability to overcome the computational difficulties arising 
when tackling large molecular systems [64J . One particular variant of subsys- 
tem DFT is the FDE approach developed by Wesolowski and Warshel [63J. 
FDE splits a system into interacting subsystems and yields subsystem elec- 
tron densities separately. Hence, it defines the total electron density, p(r), 
as a sum of subsystem densities [63] . 

# of subsystems 

p(r)= M r ). ( 8 ) 

i 

where the sum runs over all subsystems. The subsystem densities are found 
by solving subsystem-specific KS equations, known as KS equations with 
constrained electron density or KSCED [65J. In these equations, the KS 
potential vks( t ) is augmented by an embedding potential t> em &(r) which in- 
cludes the electrostatic interactions taking place between electrons and nuclei 
of the subsystems as well as a potential term deriving from the non-additive 
kinetic energy (T s nadd ) and non-additive exchange-correlation (E££ ) func- 
tional derivatives. We refer to Refs. [SHJESlEn] for more details on the FDE 
theoretical framework, however, for sake of clarity we report here the spin 
KSCED equations which lead to the subsystem orbitals 
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= 00, (9) 
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where 0£\ are the molecular orbitals of subsystem I and of spin a. 

The FDE scheme greatly reduces the computational cost compared to KS- 
DFT, as there is no need to calculate orbitals for the total ("super") system, 
and the total computational complexity is then linear over the number of 
subsystems composing the total supersystem. However, the reduced cost 
in FDE compared to KS-DFT is achieved at the expense of approximating 
the non-additive functionals with semilocal density-functional approximants. 
This approximation, for the kinetic energy especially, is the single source of 
certain shortcomings of FDE, for example when applied to covalently bound 



subsystems [67J. 

4 Charge Transfer Excitations with FDE 



Recently, two of the present authors have shown in Ref. [62] that FDE can 

be successfully used to calculate migration CT couplings and excitation ener- 
gies. In this work, we present a different algorithm for obtaining the elements 
of the full-electron Hamiltonian and overlap matrices defined in Eqs. (jlj) and 
(EJ) which has significant advantages over the previous method. Specifically, 
the new formalism can be applied to symmetric systems (out of reach be- 
fore due to a singularity in the working equations [201 [62]) an d, as it will 
be clear in Section [5j yields very accurate excitation energies with main- 
stream GGA-type functionals, while before it was necessary to include HF 
non-local exchange to the functional to counteract the self-interaction error. 
In addition, two new algorithms have been implemented here. We call them 
subsystem-joint transition density (JTD) and subsystem-disjoint transition 
density (DTD) formalisms. While the former is computationally straight- 
forward, it formally does not scale linearly with the number of subsystems. 
Instead, the latter scales linearly with the number of subsystems and, thus, 
it is theoretically consistent with the FDE formalism. 
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4.1 Approximate Couplings and Total Energies in Sub- 
system DFT 

Suppose we have two broken-symmetry Slater determinants describing two 
charge- localized states $1 and $2- The wave functions in terms of the two 
set of molecular orbitals take the form 

(10) 

where the antisymmetrizer A also contains appropriate normalization con- 
stants. Generally, the two sets of molecular orbitals {4>^} and {0^} may not 
be orthonormal to each other and within the sets. This usually happens for 
broken-symmetry HF and KS-DFT solutions as well as for KS-determinants 
derived from Constrained DFT [13] and subsystem DFT calculations [62J. 
We define the transition orbital overlap matrix as follows 

(s (12) )H = (tf ) i0! 2) ). (ii) 

In what follows, we only consider the case of a subsystem DFT calculation, 
i.e. we only deal with Slater determinants of the system constructed form 
subsystem molecular orbitals as done in a previous work [62J. To avoid 
redundance with the theory of Ref. [62], let us briefly state that, similarly 
to the Valence Bond theory [57], the subsystem DFT version of the Slater 
determinant in Eq. (fit)]) features products of occupied subsystem orbitals 
regardless of the fact that the orbitals between subsystems might not be 
orthogonal to each other. In the case of a two-subsystem partitioning, the 
S^ 12 ) transition orbital overlap matrix can be formally written as 

/ g(12) g(12) \ 

S<12 ' = { sfe s|f J ' (12) 

where and S^P are transition orbital overlap matrices calculated with 
the orbitals belonging to subsystem I or II, respectively (subsystem tran- 



$, = A 
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sition orbital overlap matrices, hereafter), while S T TI and S IZ j include the 
overlaps of the orbitals belonging to subsystem I with the ones belonging to 
subsystem II. Let us clarify that there are two sources of non-orthogonality 
in this formalism. The first one stems from the overlap between the full- 
electron charge-localized states, i.e. Eq. (JSJ), and the second one arises at the 
(subsystem) molecular orbital level in Eqs. (flOj) and (ITT]) . The orbital over- 
lap matrix is generally non-diagonal, also when it is computed from orbitals 
belonging to a single state, namely, 

(S w )« = (tf|#>. (13) 

It can be seen that the above matrix also is non-diagonal as the orbitals 
{<fi£ }k=i,N are borrowed from an FDE calculation which does not impose 
orthogonality to orbitals belonging to different subsystems [63||68]. 

Going back to the elements of the and S^ 2 ) submatrices, it is clear 
that they are small in magnitude compared to the subsystem transition over- 
lap matrices provided that the subsystems are not covalently bound to each 
other. The matrix form in Eq. (11 2p partially loses its block form in the case of 
electronic transitions that change the number of electrons in the subsystems, 
such as the ones considered in this work. Specifically, if the CT transition 
involves one electron, then the transition overlap matrix will have one col- 
umn (row) with sizable non-zero elements across the subsystems involved in 
the CT event. 

In this work, we use the following formulas for the calculation of the 
Hamiltonian coupling between the two charge-localized states [T4"tl(3§] 

H 12 = |$ 2 > = E [/9 (12) (r)] S 12 (14) 

where S\2 — det (S^ 12 ^), and 

P (12) (r) = Etf ) W( S(12) ) fc 7tf ) W (15) 

kl 
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is a scaled transition density derived from the orbitals of the two states, and 
the functional E is an appropriate density functional. The transition density 
in Eq. f|T5|) is obtained from the integration of the one-particle Dirac delta 
function, namely [70] 



($i|iV5(ri -r)|$ 2 > = N J r l •r n & x {Fi • ■ " *n)K t \ ~ r )$2(ri ■ ■ • r^) 
= ^D( 12 )(fc|/)^(r)^(r) 

= det(S( 12 ))^tf)(r)(S^) fc - 1 ^(r) (16) 

fcZ 

where D^ 12 '(k\l) is the signed minor (or cofactor) of the transition orbital 
overlap matrix of Eq. (fTTT) . Then Eq. (|T5|) is obtained by simply dividing 
Eq. (|T6|) by det (S^ 12 - 1 ), so that the scaled transition density integrates to the 
total number of electrons (N) rather than to iVdet (S^ 12 -*). 

Scaling the transition density is a matter of algebraic convenience. Eq. 
( I14p is particularly simple in terms of a the scaled transition density, and is 
rigorous if the wavefunctions are single Slater determinants and the Hamil- 
tonian is the molecular electronic Hamiltonian [6l? }ITiT[T3"] (i.e. as in the HF 
method). In our work, we consider two charge- localized states per calcula- 
tion. These are generally non-orthogonal. However, for some systems, it is 
possible that the two states could accidentally be orthogonal. In this sce- 
nario, the scaled transition density is recovered by employing the Penrose 
inverse of the transition overlap matrix in Eq. (TT5]) . The matrix is inverted 
in the N — x dimensional subspace (usually x — 1) where is it is not singular. 

It is important to point out that the coupling formula in Eq. (|14p was 
derived assuming that $1/2 are broken-symmetry HF wave functions. Even 
though the HF and KS wave functions are both single Slater determinants, 
the formula is not formally transferable from the HF method to a DFT 
method (as we do in this work). Applying Eq. (fl4"|) in the context of DFT is 
an approximation. It can be shown that in the context of linear-response TD- 
DFT Eqs. fll4p and (TO)]) are equivalent to the Tamm-Dancoff approximation 
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to TD-DFT. In addition, let us clarify that formally the HF exchange and 
HF and KS kinetic energy terms are computed with the transition density 
matrix, and thus the density in Eq. (fT5l) should be replaced by a quantity 
that depends on one electronic coordinate r as well as another electronic 
coordinate, r', namely // 12 )(r,r') = Ylu ^k\ T ) i^^) ki ^P^ 1 "')- However, 
for sake of clarity, here we will limit ourselves to reporting the (transition) 
density without introducing the density matrix notation. As an example 
of how we apply Eq. (fllj) in practical calulations, consider a calculation in 
which we employ the LDA exchange density functional, then the exchange 
contribution to the off-diagonal Hamiltonian matrix element becomes 

E x [p^(r)j S 12 = -Sj- (i)' J (p< 12 >(r))*r . (17) 

The diagonal elements of the Hamiltonian are computed in a similar fash- 
ion as the off-diagonal ones, namely 

H H = ($ i \H\Q i )=E[p®(T)] (18) 

where 

P W (r)=E^( r )( S(i °)« ( 19 ) 

kl 

Note that the overlap element has disappeared in Eq. ([TBI . While this is a 
trivial consequence of the normalization condition on the Slater determinant 
in a regular KS-DFT calculation, the KS Slater determinant of the super- 
system built with subsystem orbitals is not normalized. However, the cor- 
responding FDE density theoretically is the correctly normalized correlated 
density of the full system, integrating to the total number of electrons |62j . 
This is an important point as the number of electrons in a calculation is set 
by the trace of the (transition or diagonal) orbital overlap matrix, while the 
norm of a Slater determinant depends upon the determinant of the corre- 
sponding orbital overlap matrix. 

To summarize, the approximations we employ are equivalent to (1) re- 
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placing the HF Slater determinants needed to define the charge-localized 
states $1/2 with Slater determinants made up of subsystem KS molecular or- 
bitals, and (2) replacing the HF expression of E [p( 12 )(r)] by an approximate 
density functional. 

Approximation (1) is often invoked in DFT calculations of electronic cou- 
plings HHT51E31ES1E2], while (2) was used before by Kaduk et al. [H], and it 
is introduced here for the first time in the context of subsystem DFT. Con- 
trary to the density, the scaled transition density introduced in Eq. (fl5|) might 
feature regions of space where it carries a negative sign. Generally, exchange- 
correlation and kinetic-energy functionals may feature arbitrary fractional 
powers of the density that would render the exchange- correlation energy 
complex. By working with a scaled transition density that is quasi-positive 
almost everywhere in space (see below), the simplest pragmatic workaround 
is to set the scaled transition density to zero everywhere in space where it 
actually has negative values. While this may seem to be a somewhat drastic 
approximation, we note that in practice only a small part of the scaled tran- 
sition density needs to be changed, and the very good results observed in our 
test cases below empirically justify this procedure (see Section [5]). In all the 
cases considered here, the neglect of the negative parts of the scaled tran- 
sition density affected its total integral by less than a tenth of a percentile 
point. 

The quasi positivity of the scaled transition density can be explained with 
two arguments. In a frozen core approximation, if the only orbital product 
term left in the scaled transition density is negative almost everywhere in 
space, then that product will give rise to a negative element of the inverse of 
the transition overlap matrix (being just the inverse of that orbital overlap 
in this approximation). Thus, the scaled transition density will be positive 
almost everywhere, as the above mentioned orbital product is multiplied by 
the inverse overlap in the equation for the scaled transition density. The 
second argument uses the fact that the scaled transition density integrates 
to N, it can be decomposed into two contributions, one corresponding to the 
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N — 1 electron system undergoing little changes in the transition, and the 
other corresponding to the single transferring electron. Only the latter com- 
ponent has possibly negative contributions to the scaled transition density, 
and is unlikely to overpower the much larger N — 1 component. 

Generally, in those cases where the S12 overlap is zero because x orbital 
overlaps (x rows/columns of the transition overlap matrix) are identically 
zero, we employ the Penrose inversion of the transition overlap matrix and 
p( 12 ) is indeed a N — x "positive" electron system "plus" an x electron system 
that overall integrates to zero (and thus half positive and half negative). 
In this case the scaled transition density integrates to N — x and because 
x <C N, the negative parts of the x electron system will unlikely overpower 
the positive N — x electron system. In any even, this orthogonal case will 
yield a zero Hyi matrix element as calculated by Eq. (114)) . Thus, even in those 
cases where the scaled transition density has large negative parts (which we 
never encountered in our calculations), its contribution to the Hamiltonian 
matrix element will be negligible due to either a negligible or a zero overlap 
element. 

We should warn the reader that the sign of the scaled transition density 
is independent of wave function phase factors. In fact, wave function phase 
factors do not affect the scaled transition density as the overlap element be- 
tween the iV-electron wave functions (Slater determinants in the context of 
our work) will carry with it the same phase factor as the conventional tran- 
sition density. Thus, following the definition of the scaled transition density 
in Eq. ( TT5l) . dividing the conventional transition density by the overlap has 
the effect of removing the phase factor altogether. 

For sake of brevity, hereafter we will drop "scaled" when referring to 
the above scaled transition density. In calculating the transition density 
from FDE subsystem orbitals coming from the KSCED equation in Eq. (JU]) 
by applying Eq. ffl5l) . one must consider all of its subsystem contributions 
as well as its overlap-mediated inter-subsystem couplings. The transition 
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density for a system partitioned in Ns subsystems becomes 

N S 

i,j kei leJ 

It is interesting to notice that the above transition density "joins" orbital 
transitions of one subsystem with the ones of another subsystem through the 
terms in the sum over the subsystem's labels where I ^ J. The coupling 
of these inter-subsystem transitions are weigthed by the matrix elements of 
the inverse of the transition orbital overlap matrix, (S^ 12 -*) . The diagonal 
density, i.e. the density of the broken-symmetry charge-localized states, can 
be formulated similarly applying Eq. f fl9|) . namely, 

N S _ 

^w-EEE^wW-CW' < 21 > 

I,J k£l l£j 

As pointed out previously [62] this density is not exactly the sum of subsystem 
densities obtained with the FDE calculation from Eq. (jSj), as the intersubsys- 
tem overlap elements do not appear in the FDE formalism. Such a density 
mismatch constitutes a second-order effect [62] and is not a concern here as 
we are going to apply the method only to non-covalently bound subsystems, 
that is, super-systems with small inter-subsystems orbital overlap. As we 
will see in the next section, this density issue is completely by-passed by an 
FDE-compatible theory we call "disjoint transition density formulation". 

4.2 Subsystem-Joint and Subsystem-Disjoint Transi- 
tion Density Formulations 

The transition density in Eq. (1201) can be approximated by a density com- 
posed only of intra-subsystem orbital transitions. By recognizing that the 
off-diagonal blocks S\ TI in Eq. ffl2l are always very small in magnitude when 
the subsystems are not covalently bound and if there is no electron transfer 
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between subsystems the disjoint transition density can be written as 

/ffi>O0 = £p? 2) 0O, (22) 
i 

k,l&I 

where the tilde has been included to distinguish the above definition of dis- 
joint transition density with the one we use for charge-transfer transition, see 
below. The off-diagonal Hamiltonian matrix element (.Hi 2) can therefore be 
calculated with either the joint transition density (JTD) in Eq. f[2"Uj) . giving 
rise to the JTD Hamiltonian couplings; or with the disjoint transition density 
(DTD) of Eq. (f22l) . giving rise to the DTD Hamiltonian couplings. 

The above definition of disjoint transition density can be cast in terms 
of a pure subsystem DFT formulation. Because in the FDE formalism both 
the ground and excited state densities can be written as a sum of subsystem 
densities, as in Eq. (JBJ), namely, 

5>P(r), and (24) 
1 

E^w- ( 25 ) 
1 

Given that the subsystem densities are derived from Slater determinants of 
subsystem orbitals, then it is enough to invoke approximation (1) (i.e. replace 
the HF Slater determinants needed to define the charge- localized states $i/ 2 
with Slater determinants made up of subsystem KS molecular orbitals) the 
transition density relating states 1 and 2 is of the DTD type [as in Eq. fT22l ], 
However, if the electronic transition involves electron transfer between sub- 
system K and L, then the transition density cannot be completely disjoint, 
and it can be approximated by (note, compared to Eq. ( 1221) the tilde is re- 



PdtdM = 
Pdtd ( r ) = 
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moved) 



Np 

/>88>M = *£?M + £ 



Pi ( r ) 



k.lEK.L 



(26) 
(27) 



where Sj 12 ^ are the subsystem transition orbital overlap matrices and is 
the full transition orbital overlap matrix of the combined K and L subsystems 
and pj{r) is the same as in Eq. fT23|) . 

The diagonal elements of the Hamiltonian over the broken-symmetry 
states in the DTD calculations become 



Hu = (H> i \H\$ i ) = E 



PdtdM 



nadd 



Pl,Pli,--- 



+ E. 



nadd 



(*) « 
Pl,PlI,--- 



(28) 



For the above expressions, we use the FDE total energy partitioning in sub- 
system energies and non-additive functionals [71]. Similarly, the off-diagonal 
elements become 



#12 = <$l|#|$2> = E 



( 12 ) / > 

PdtdI 1- ) 



5 



5* 



12 



12 



nadd 



(12) (12) 
Pj >P// . 



nadd 



(12) (12) 

Pi ,Pii > 



(29) 



where the sums run over the number of disjoint subsystems in the calculation 
and the non-additive functionals have the same meaning as in the regular 
FDE formalism [5U[63l[65],[66] ■ Note that for charge-transfer transitions, the 
above sum over the disjoint subsystems is carried out in a similar way to 
the one in Eq. ( 126]) . i.e. the subsystems undergoing variation of the number 
of electrons [K and L in Eq. (J27J] are grouped together in one supersystem 
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that maintains an overall constant number of electrons during the transfer 
of charge. 

Let us remark that the computational cost needed to calculate the matrix 
elements in the DTD formalism is linear- scaling with the number of subsys- 
tems, while the JTD formalism is not - it scales as N 3 with iV being the 
number of electrons in the supersystem even though with a much smaller 
scaling coefficient than regular KS calculations (as in the JTD formalism the 
KS equations do not need to be self-consistently solved for the supermolec- 
ular system. Instead the transition orbital overlap matrix is inverted only 
once) . 

5 Validation of the methodology 

After implementation of the above described JTD and DTD algorithms 
in ADF, in what follows, we aim at providing a thorough benchmark of 
the method. We compare our calculated migration CT excitation energies 
against two dimer systems: a water dimer radical cation, (H^O)^, and an 
ethylene dimer radical cation, (C2H4)^~, at varying inter-monomer separations 
for which we carried out several different high-level wave function method cal- 
culations of the excitation energies. We also calculate hole transfer excitation 
energies for several DNA nucleobase dimers, and compare them to CASPT2 
calculations. The idea behind this benchmark is to seek a validation of the 
ability of our FDE method (note JTD and DTD are equivalent for only two 
subsystems) to reproduce quantitatively (with at least chemical accuracy, 1 
kcal/mol or 0.04 eV deviation) the CT excitation energies calculated with 
the wave function methods considered. Unless otherwise noted, all DFT 
calculations have been carried out with the PW91 density functional and a 
TZP (triple-zeta with polarization) basis set with ADF, all EOM-CCSD(T) 
calculations have been carried out with the NWChem 5.2 software [75], and 
all Fock-Space CCSD calculations have been carried out with the DIRAC 
11 software [76J. 
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5.1 Water Dimer 

The nuclear positions of a neutral water dimer were first optimized with a 
DFT method employing the BP86 functional and a TZP basis set (triple 
zeta with polarization). Figure SI [77] depicts the obtained structure. The 
subsequent calculations were carried out on the same system by translating 
one water molecule away from the other along the hydrogen-bond axis. The 
results are summarized in Figure [TJ 

[Figure 1 about here.] 

We calculated the excitation energy in two ways. In the first one (FDE/FDE 
label in the figure) the diabatic energies [or the diagonal elements of the 
diabatic Hamiltonian in Eq. (j3J)] were computed using Eq. ( 128]) with the 
FDE density calculated as in Eq. ([8]). In the second calculation (FDE/ET) 
the diabatic energies were calculated by integrating the energy functional as 
in Eq. (lisp with the density computed from the FDE subsystem orbitals as 
in Eq. (TT9]) . 

We compare our values with excitation energies from four coupled cluster 
calculations. Three of these calculations were carried out with the EOM-CC 
method (with different basis sets and reference determinant). In the fourth 
calculation the excitation energy is obtained by subtracting from each other 
the ionization potentials obtained with Fock-Space Coupled-Cluster [78J [FS- 
CCSD] calculated from two different guesses of the initial ionized state: one 
being the HOMO, almost entirely localized onto one monomer, and the other 
one the HOMO-1, almost entirely localized onto the other monomer. 

The comparison of the CC methods with our values is very good for the 
FS-CCSD values while the EOM-CCSD(T) values are consistently larger in 
magnitude. This behavior of EOM-CC methods has been noted before for 
MCT excitations [19] and, as mentioned in the introdution, it is due to the 
fact that the EOM-CC excited-state ansatz cannot be written in the regular 
exponential form, thus making EOM-CC non-size extensive jl3]. One can 
rationalize this with an argument based on orbital relaxation. In the excited 
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state, the charge has migrated to another location and the orbital of the 
donor must relax back to the neutral state. Similarly, the orbitals of the 
acceptor must relax to the charged state. In both cases single excitations 
are the main contribution, however they must be computed on top of the 
single excitation needed to describe the CT itself. Therefore one needs at 
least double excitations to describe these relaxation effects [22] • Adding 
perturbatively the triple excitations usually further improves this situation 
allowing the EOM-CCSD(T) to reach typical deviations of 0.1-0.3 eV on 
excitation energies [49J. 

We therefore consider the EOM-CCSD(T) energy overestimation a good 
sign that our method yields accurate excitation energies. Convincing evi- 
dence is provided by the FS-CCSD energies, as they compare with our values 
to within chemical accuracy at every inter-monomer separation. An issue of 
orthogonality arises, as the ionization potentials calculated by the FS-CCSD 
method are deduced from cationic wave function which are not strictly or- 
thogonal to each other. This however, plays a very minor role here, as the 
overlap and the electronic coupling between the two diabatic states is al- 
most zero as proved by our calculations of S12 in Table SI and S2 in the 
supplementary materials [77]. 

For sake of completeness, in the supplementary materials [77] we report 
the numerical values used to obtain the plot in Figure [1] as well as additional 
FDE calculations carried out with the BLYP functional showing a similar 
behavior to the PW91 calculations. 

5.2 Ethylene Dimer 

The nuclear positions of a neutral ethylene molecule were first optimized 
with a DFT method employing the BP86 functional and a TZP basis set 
(triple zeta with polarization). Then a copy of the same molecule is pasted 
so as to obtain a 7r-stacked ethylene dimer with intermolecular separation of 
R. Figure S2 [77] depicts the obtained structure. The calculated excitation 
energies are plotted in Figure EJ 
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[Figure 2 about here.] 

Due to the symmetric arrangment of the ethylene dimer, this system be- 
haves in a very different way from the water dimer considered above. Here, 
the orbital relaxation issues faced before are less important as the hole is 
completely delocalized between the two monomers - the inner orbitals are 
not dramatically affected by the excitation. The hole sits in orbitals with 
different symmetry in the ground and in the CT excited state. The ground 
state has a symmetric HOMO having a g symmetry, while the first excited 
state has the HOMO of b 2u symmetry. Because of these reasons, we expect 
the EOM-CC calculations to be of much higher accuracy than for the wa- 
ter dimer case. And in fact this is what we notice. Our calculations are in 
very good agreement with both the EOM-CCSD(T) and the FS-CCSD cal- 
culations. Our values slightly underestimate the CC calculation at shorter 
separations, even if only of 0.1 eV. At larger separations (R> 6), due to the 
lowering magnitude of the excitation energy, and to a much reduced extent 
compared to the water dimer case, the overestimation of the excitation en- 
ergy by EOM-CCSD(T) becomes noticeable again. Instead, FS-CCSD and 
our calculations are still in very good agreement. The EOM-CC calculations 
carried out with the 6-311G(d,p) basis set show a slight off-trend behavior. 
We attribute this to the absence of augmented functions in the basis set. 

Once again, for sake of completeness, in the supplementary materials [77] 
we report the numerical values used to obtain the plot in Figure [2] as well as 
additional FDE calculations carried out with the BLYP functional showing 
a similar behavior to the PW91 calculations. 

5.3 Hole transfer in DNA tt— stacked nucleobases 

In a recent paper [62] we have applied FDE to constrain charges on DNA 
nucleobases and subsequently we calculated electronic couplings and CT ex- 
citation energies. In the previous work we employed the formalism of Migliore 
et al. [2"m i53 | IT9] for the calculation of the coupling which retained all the de- 
ficiencies of the approximate density functionals, particularly self interation. 
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The self interaction error caused a gross overestimation of the coupling and 
excitation energies which was ameliorated by mixing in non-local Hartree- 
Fock exchange in the functional. 

Here we recalculated the same systems, i.e. Guanine dimer, (GG) + , and 
the Guanine-Thymine complex, (GT) + . In addition, here we report cal- 
culations on the Adenine dimer, (AA) + , Thymine dimer, (TT) + , Thymine- 
Guanine, (TG) + , Adenine- Guanine, (AG) + , and the Guanine-Adenine, (GA)" 1 
The results are summarized in Table [TJ 

[Table 1 about here.] 

Comparing the CT excitation energies calculated with the method of Migliore 
with the ones calculated in this work, it is clear that the self interaction error 
does not affect the new calculations. This can be understood by consider- 
ing that the formalism of Migliore requires the KS wave function of the 
supermolecular system (i.e. composed of the donor and acceptor). The self 
interaction error, even though it is present in all calculations, does not affect 
the FDE calculations because the hole there is localized by construction and 
is not allowed to overdelocalize [62|l66]. 

The excitation energy values are generally in good agreement with the 
CASPT2 benchmark values. For the (TT) + system no high-level bench- 
mark calculations are available, while for the (AA) + and (TG) + systems 
our excitation energies deviate by about 0.1 eV against the CASPT2 val- 
ues. After inspecting the excitation energies calculated with the progres- 
sion CASSCF(7,8), CASSCF(11,12) and CASPT2(11,12) in Ref. [SO] for the 
(AA) + and (TG) + systems we notice that while for all the other nucleobase 
stacks the excitation energies calculated with the three methods are similar 
to each other, the cases of (AA) + and (TG) + stand out. For example the 
CASSCF(7,8) excitation energies are 0.144 eV and 1.235 eV for the (AA)+ 
and (TG) + systems respectively. Then these values drop to 0.047 eV and 
1.097 eV for the CASSCF(11,12) and then back up to 0.097 eV for the (AA) + 
system but go down to 0.797 eV for (TG) + in the CASPT2 calculation. 
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According to Blancafort and Voityuk [80J, these large fluctuations in the 
excitation energy are to be expected especially when going from CASSCF to 
CASPT2. However, they also notice that for the (AA) + system a larger active 
space than they employed should be considered to confirm the couplings and 
excitation energies for this nucleobase dimer complex. Therefore, it is difficult 
to compare the excitation energies for all the nucleobase stacks we compute 
here with the CASSCF and CASPT2 calculations of Blancafort and Voityuk. 
Further analysis of these calculations (i.e. study the effect of the size of the 
active space in the CASSCF/CASPT2) lie beyond the scope of this work. 

6 Pilot calculations: embedding effects in wa- 
ter and ethylene clusters 

In this section we present calculations carried out for systems composed of 
more than 2 subsystems. As pilot studies we choose the same small dimer 
systems considered above and embed them in clusters of the same molecule. 
First we consider the hole transfer in ethylene dimer embedded in an ethylene 
matrix, then we consider a water dimer in a cluster extracted from liquid 
water. 

6.1 Ethylene clusters 

In Figure [3] we depict a selected number of ethylene clusters used in the 
calculation, and in Table [2] we report the values of excitation energy and 
electronic coupling corresponding to the hole transfer in the selected dimer. 
The dimer has the same geometry of the dimer considered in Section I5T21 with 
inter-monomer separation of R = 4.0 A. For the DTD calculation the non- 
additive exchange-correlation functional used is PW91, while PW91k [8T||82] 
was employed for the evaluation of the non-additive kinetic-energy functional. 

[Figure 3 about here.] 



25 



[Table 2 about here.] 



The excitation energy and coupling values reported in Tableware essentially 
independent of the size of the cluster, showing that for this dimer system 
we should expect a minor environmental effect of the CT excitation and 
the hole transfer electronic coupling. This is indeed what was found for 
this conformation and inter-monomer orientation in a study by Lipparini 
and Mennucci |83j where the environment was modelled as a polarizable 
dielectric. 

It is important to notice how the JTD and DTD calculations yield very 
similar excitation energies and couplings. The DTD values always lie a few 
meV lower than the JTD ones. The RMS deviation of the excitation energies 
calculated with the JTD versus DTD method is 1.5 meV, while the average 
deviation is 1.2 meV. While it is difficult to pinpoint the origin of this small 
discrepancy, we note that the contribution to the excitation energy and elec- 
tronic coupling from the non-additive kinetic energy functional ranged from 
only a few meVs in most cases to 23 meV in the four-member cluster system. 
In all cases, this contribution brought the excitation energy value closer to 
the JTD one. 

Lipparini and Mennucci [53], also considered dimers with different orien- 
tation by rotating one of the ethylene molecules around the carbon-carbon 
axis. In the supplementary materials [77], we report calculations carried out 
for the dimer systems with a rotated monomer in the presence of a minimal 
environment. Our calculations largely reproduce the findings of Lipparini 
and Mennucci. However, our FDE calculations allow for the recovering of 
effects of the discrete molecular environment. 

6.2 Water clusters 

The clusters were generated by chosing a water dimer complex from bulk wa- 
ter and then carving the desired water molecules around it. The coordinates 
of a large cluster of bulk water was given to us by Daniel Spangberg [84|. 
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More details on the generation of the bulk water is available in the supple- 
mentary material [77]. 

Figure 0] depicts a selected set of clusters considered with the highlighted 
water dimer chosen for the hole transfer. All the other molecules belong to the 
environment and polarize according to where the hole is located. However, 
they do not undergo variation of electron number during the hole transfer 
process. In Table |3]we collect the values of electronic couplings and excitation 
energies calculated with the JTD and DTD formalisms. As noted in the 
ethylene dimer calculations, also here the DTD values always underestimate 
the JTD ones. In this case by 1-60 meV for excitation energies and 1-15 meV 
for the electronic couplings. The RMS deviation of the excitation energies 
calculated with the JTD versus DTD method is 0.06 eV, while the average 
deviation is 0.05 eV. 

In our water cluster calculation we see a very different behaviour from the 
ethylene cluster. The electronic couplings are as effected by the environment 
as the excitation energies. This is likely the effect of water having a perma- 
nent dipole, thus adding large contributions to the electric field interacting 
with the hole. Contrary, in the ethylene clusters case the environmental elec- 
tric field acting on the hole is of much weaker magnitude as it is due only to 
polarization of the ethylene molecules. What we notice in our calculations is 
that the more water molecules are included in the cluster, the more the elec- 
tronic coupling decreases. This was also a feature of the effect of embedding 
by water on excitonic couplings in 7r-stacked chromophore diads [85] . For the 
configuration considered, the excitation energy decreases by about 1 eV as 
the cluster size increases. Such large variations of the excitation energy in 
condensed phases with polar solvents is well known [86] . Witness of this are 
also the much larger reorganization energies of electron transfer processes in 
polar solvent than in non-polar ones [3"| l8"71 - [50"] . 

[Figure 4 about here.] 

[Table 3 about here.] 
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7 Conclusions 



In this work, we have developed a new and linear-scaling DFT method aimed 
at accurately predicting charge-transfer excitation energies and diabatic cou- 
plings among charge-localized states, and validated it against high-level wave 
function calculations. 

The success of the disjoint transition density method (or DTD, see sec- 
tion U]2]) resides in two important properties. First, it computationally scales 
linearly with the number of non-covalently bound molecules in the system. 
Secondly, the hole transfer excitation energies calculated with mainstream 
GGA-type functionals reproduce within chemical accuracy calculations car- 
ried out with several types of Coupled Cluster methods. For the benchmark 
cases considered, our method outperforms the biased EOM-CCSD(T), while 
it is in excellent agreement with vertical ionization potential differences cal- 
culated with Fock-Space CCSD(T) calculations. We also carried out calcu- 
lations of the hole transfer excitations in seven DNA base pair combinations 
reproducing CASPT2(11,12) calculations within a 0.1 eV deviation or bet- 
ter. Pilot calculations on molecular clusters of water and ethylene have also 
been carried out with cluster sizes up to 56 molecules. We show how the 
method is able to recover solvation effects on diabatic couplings and site 
energies (needed in electron transfer calculations) as well as the excitation 
energies at the discrete molecular level and fully quantum-mechanically. This 
constitutes an important step forward in the development of linear scaling 
electronic structure methods tailored to electron transfer phenomena. 

Future improvements of the method are underway and will include the 
implementation of hybrid functionals in the post-SCF calculation of the full- 
electron Hamiltonian matrix elements and the extension of the method to 
both excitations involving charge separation and covalently bound subsys- 
tems. 
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Figure 1: Migration charge-transfer excitations (in eV) for water dimer radi- 
cal cation in the C s symmetric configuration with an intermonomer displace- 
ment from the equilibrium distance of < R < 15 A. 
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Figure 2: Migration charge-transfer excitations (in eV) for ethylene dimer 
radical cation in the D 2 h symmetric configuration with an intermonomer 
displacement from the equilibrium distance of 6 < R < 15 A. Full 3 < R < 15 
A range available in Table S3 [77]. 
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Figure 3: Depiction of the ethylene clusters used in the calculations contain- 
ing: inset (a) 4, (b) 6, (c) 8, and (d) 10 ethylene molecules. The depicted 
clusters were obtained from the 20-molecule cluster cutting the furthest away 
molecules from the center of mass of the two molecules experiencing the hole 
transfer. Figure obtained with molden |91j . 
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Figure 4: Depiction of the water clusters used in the calculations containing: 
inset (a) 2, (b) 4, (c) 6, (d) 8, and (e) 10 water molecules. The depicted 
clusters were obtained from the 56-molecule cluster cutting the furthest away 
molecules from the center of mass of the two molecules experiencing the hole 
transfer. Figure obtained with vmd [H2] and molden 
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Table 1: Hole transfer excitations and couplings in 7r-stacked DNA nucle- 
obase radical cation complexes. All energy values in eV. The - stands for 
"not available". 



Composition 


Sl2 


V 12 










(AA)+ 


0.004 


0.092 


0.198 






0.097 


(AG)+ 


0.002 


0.177 


0.421 






0.340 


(GA)+ 


0.012 


0.058 


0.530 






0.560 


(GG)+ 


0.009 


0.051 


0.405 


1.653 


0.418 


0.392 


(GT)+ 


0.020 


0.104 


1.082 


1.944 


1.159 


1.175 


(TG)+ 


0.009 


0.056 


0.657 






0.797 


(TT)+ 


0.018 


0.099 


0.208 









a This work. 

b PW91/TZP from Ref. [62] 
c BHandH/TZP from Ref. [62] 
d CASPT2 values form Ref. [SD] 
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Table 2: Ethylene cluster electronic couplings and charge transfer excitations. 
The hole transfer occurs between two ethylene molecules separated by R = 
4.0 A. All values in eV. N is the number of molecules in the cluster. 



JTD DTD 



N 


V 12 


A£ ex 


v 12 


A£ cx 


2 


0.260 


0.521 






4 


0.261 


0.539 


0.261 


0.540 


6 


0.262 


0.524 


0.260 


0.521 


8 


0.262 


0.535 


0.261 


0.534 


10 


0.262 


0.538 


0.260 


0.538 


20 


0.262 


0.535 


0.260 


0.534 
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Table 3: Water cluster electronic couplings and charge transfer excitations. 
All values in eV. N is the number of molecules in the cluster. 



JTD DTD 

N V 12 AE CX V 12 AE ( 



Z 


n 
u 


.zoy 


i 
i 


.4ol 










3 


0. 


.242 


1 


.138 





.253 


1 


.098 


4 





.250 


1 


.447 





.249 


1 


.422 


5 





.226 


1 


.254 





.241 


1 


.189 


6 





.251 





.951 





.250 


0, 


.913 


7 





.227 





.771 





.242 


0, 


.716 


8 





.257 


1 


.527 





.252 


1 


.498 


9 





.236 


1 


.344 





.244 


1 


.275 


10 


0. 


.261 


1 


.416 


0. 


,256 


1 


.383 


11 


0. 


.248 


1 


.252 


0. 


,254 


1 


.187 


13 


0. 


.249 


1 


.227 


0. 


,255 


1 


.154 


17 


0. 


.245 


1 


,094 


0, 


.251 


1 


.013 


21 


0. 


.211 


1 


.002 


0, 


,219 


0, 


.929 


27 


0. 


.160 





.610 


0, 


.170 


0, 


.544 


31 


0. 


.177 





,743 


0, 


,186 


0, 


.680 


36 


0. 


196 





,843 


0. 


.203 


0, 


.777 


41 


0. 


.199 


0. 


.733 


0. 


,205 





.700 


46 


0. 


.155 





.563 


0. 


,165 





.518 


51 


0. 


.147 


0. 


.607 


0. 


,160 





.558 


56 


0. 


.128 


0. 


.468 


0. 


,143 





.422 
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